      PROGRAM NAD3(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,PUNCH)
C                            PROGRAM NAD3
C
C NAD3 IS A FINITE DIFFERENCE PROGRAM EMPLOYING MURRAY LANDIS TRANSFORMATION
C FOR THREE PHASE BINARY SYSTEMS. THIS PROGRAM ALLOWS FOR PLANAR, CYLINDRICAL,
C OR SPHERICAL GEOMETRY. B STANDS FOR THE SOLUTE OR FILAMENT ELEMENT WHILE A
C STANDS FOR MATRIX ELEMENT. ALL CONCENTRATIONS ARE READ IN ATOM FRACTION B,
C ALL DISTANCES IN CENTIMETERS AND TIME IN SECONDS.
C
      REAL KPRIME
      INTEGER PCS
      COMMON/STORE/NX,CN( 51),DN( 51),ND,CP(501),DP(501),DMAX
      DIMENSION XCAP(7),YCAP(7),TYTLE(8),RESULT(3),XN(501)
      DIMENSION FI( 8),ABLE(16),XSN(501),CSN(501)
      DATA RESULT/10HFID FOR TH,10HREE PHASE ,10HDIFFUSION /
C
C PSEUDO AND FONTS ARE SYSTEM SUBROUTINES NECESSARY FOR PLOTTING.
C
      CALL PSEUDO
      CALL FONTS(1)
C
C                             FORMAT STATEMENTS
C
  100 FORMAT(2X,//)
  101 FORMAT (3I4)
  102 FORMAT(8A10)
  103 FORMAT(5E14.5)
  104 FORMAT ( 1H1,//)
  105 FORMAT ( 1H , 8A10)
  106 FORMAT(2X,10I8)
  107 FORMAT(//2X, 36HGRID CHANGE WAS INITIATED WHEN TIME= ,E11.4,9H SEC
     1ONDS. /2X,55HFOLLOWING DATA PERTAINS TO SOLUTION BEFORE GRID CHANG
     2E. /2X,22HBETA-PHASE THICKNESS = ,E11.4,13H CENTIMETERS./2X,
     3 22HGAMMA-PHASE THICKNESS= ,E11.4, 13H CENTIMETERS./2X,
     422HALPHA-PHASE THICKNESS= ,E11.4,13H CENTIMETERS. /)
  108 FORMAT(//2X,E12.4,5X,16HATOMIC MASS OF A/2X,E12.4,5X,16HATOMIC MAS
     1S OF B /2X,E12.4,5X,20HPAR. MOLAR VOL. OF A /2X,E12.4,5X, 20HPAR.
     2MOLAR VOL. OF B //)
  110 FORMAT(F8.5,2X,E14.8)
  111 FORMAT(2X,I4,2E15.6,2X,I4,2E15.6,2X,I4,2E15.6)
  113 FORMAT(4X,13HCONCENTRATION ,5X,21HDIFFUSION COEFFICIENT/(3X,F11.4,
     112X,E12.4))
  114 FORMAT(//5X,I4,13X,72HNXB-NUMBER OF X VALUES TO BE SET UP IN THE B
     1ETA PHASE OF DIFFUSION ZONE /5X,I4,13X, 72HNXG-NUMBER OF X VALUES
     2TO BE SET UP IN THE GAMMA PHASE OF DIFFUSION ZONE/5X, I4,13X, 72HN
     3XA-NUMBER OF X VALUES TO BE SET UP IN THE ALPHA PHASE OF DIFFUSION
     4 ZONE /2X,E15.6,5X,
     527HTHB-THICKNESS OF BETA PHASE/ 2X,E15.6,5X, 28HTHG-THICKNESS OF G
     6AMMA PHASE/2X,E15.6,5X, 28HTHA-THICKNESS OF ALPHA PHASE
     7 /2X,E15.6,5X, 25HTOTAL LENGTH OF DIF. ZONE/2X,E15.6,5X,
     830HFINAL LENGTH OF DIFFUSION ZONE /2X,E15.6,5X,50HINITIAL ESTIMATE
     9D DIFFUSION DISTANCE IN BETA PHASE /5X,F4.1,13X,3HPCS)
  115 FORMAT ( //3( 3X, 3HNO., 5X, 8HDISTANCE,  5X, 11HCOMPOSITION,1X))
  116 FORMAT(//1H ,5X,23HSTARTING DIFFUSION TIME , 7X,1H=,E15.8,1X,
     16HSECS =  ,E15.8, 5H MINS /6X,20HFINAL DIFFUSION TIME,
     210X, 1H=, E15.8,    7H SECS = ,E15.8, 5H MINS
     3/6X,18HINITIAL THICKNESS=, E11.4,17H FOR BETA-PHASE, ,
     4E11.4,18H FOR GAMMA-PHASE, , E11.4,16H FOR ALPHA-PHASE
     5/6X,18HFINAL THICKNESS  =, E11.4,17H FOR BETA-PHASE, ,
     6E11.4,18H FOR GAMMA-PHASE, , E11.4,16H FOR ALPHA-PHASE
     7/6X, 24HLENGTH OF DIFFUSION ZONE ,1H=,E15.8,2X/ 6X,
     842HTOTAL DIST. BETA/GAMMA INTERFACE HAS MOVED, 1X,1H=,E15.8 / 6X,
     A43HTOTAL DIST. GAMMA/ALPHA INTERFACE HAS MOVED,    1H=,E15.8 /
     B6X, 25HINITIAL AREA UNDER CURVE= ,E15.8/6X, 25HFINAL AREA UNDER CU
     CRVE  = ,E15.8 //)
  118 FORMAT(//
     45X,F10.5,10X,46HCBG-SOL. LIMIT IN BETA AT BETA/GAMMA INTERFACE/
     55X,E13.5,5X,23HDBG-DIF. COEFF. FOR CBG //
     65X,F10.5,10X,47HCGB-SOL. LIMIT IN GAMMA AT BETA/GAMMA INTERFACE /
     75X,E13.5,5X,23HDGB-DIF. COEFF. FOR CGB //
     85X,F10.5,10X,48HCGA-SOL. LIMIT IN GAMMA AT GAMMA/ALPHA INTERFACE /
     95X,E13.5,5X,23HDGA-DIF. COEFF. FOR CGA //
     A5X,F10.5,10X,48HCAG-SOL. LIMIT IN ALPHA AT ALPHA/GAMMA INTERFACE /
     B5X,E13.5,5X,23HDAG-DIF. COEFF. FOR CAG ,2(/))
  120 FORMAT (3X,I3,10X,E15.8,10X,E15.8)
  121 FORMAT(10F6.2)
  122 FORMAT(I2,2X,7A10)
  124 FORMAT (6F10.5)
  125 FORMAT(4E15.6)
  132 FORMAT(//2X,27HINITIAL COMPOSITION PROFILE )
  133 FORMAT(//2X,25HPROFILE AFTER GRID CHANGE )
C
C            PLOT PARAMETERS
C
C XAL    LENGTH OF X-AXIS IN INCHES.
C YAL    LENGTH OF Y-AXIS IN INCHES.
C DXN    NUMBER OF INTERVALS ON X-AXIS WHERE THE TICK MARKS ARE TO BE PLACED
C DYN    NUMBER OF INTERVALS ON Y-AXIS WHERE THE TICK MARKS ARE TO BE PLACED
C XMIN   VALUE CORRESPONDING TO FIRST TICK MARK ON THE X-AXIS
C YMIN   VALUE CORRESPONDING TO FIRST TICK MARK ON THE Y-AXIS
C XIN    SPACING BETWEEN SUCCESSIVE TICK MARKS ON THE X-AXIS
C YIN    SPACING BETWEEN SUCCESSIVE TICK MARKS ON THE Y-AXIS
C XALM   MAX LENGTH OF X-AXIS IN INCHES.
C SPS    OPTION FOR SAVING PLOT SPECIFICATIONS.
C        0, AUTO SCALING BY INCREASING XIN SO THAT THE MAXIMUM
C           DISTANCE XMAX IN THE PROFILE IS LESS THAN XALM.
C        1, SAVE THE PLOT SPECIFICATIONS BY IGNORING THE POINTS
C           FOR WHICH X IS GREATER THAN THE LAST X-TICK MARK.
C NCX    NUMBER OF CHARACTERS IN THE X-CAPTION.
C XCAP   X-CAPTION
C NCY    NUMBER OF CHARACTERS IN THE Y-CAPTION.
C YCAP   Y-CAPTION
C TYTLE  TYTLE TO BE PRINTED IN THE PLOT.
C NA     PARAMETER DICTATING NUMBER OF SOLUTIONS PLOTTED ON EACH GRAPH.
C        1, ONE SOLUTION PER GRAPH
C        2, PLOT ALL SOLUTIONS ON THE SAME GRAPH.
C
      READ(5,121) XAL,YAL,DXN,DYN,XMIN,YMIN,XIN,YIN,XALM,SPS
      READ(5,122)NCX,XCAP
      READ(5,122) NCY,YCAP
      READ(5,102) TYTLE
      READ(5,101) NA
C
C ABLE      IDENTIFICATION FOR SOLUTION.
C PCS       PARAMETER DICTATING FORM OF SOLUTION
C           0 , PLANAR GEOMETRY
C           1 , CYLINDRICAL GEOMETRY
C           2 , SPHERICAL GEOMETRY
C NSP       OPTION FOR READING STARTING CONCENTRATION PROFILE
C           1, READ AN INITIAL PROFILE
C           0, THE PROGRAM INITILISES THE PROFILE WITH CPRIME IN THE
C           FILAMENT, CZERO IN MATRIX, AND A LINEAR PROFILE IN GAMMA.
C NO        PARAMETER DICTATING FORM OF OUTPUT.
C           1, WRITTEN ONLY
C           2, WRITTEN AND PUNCHED
C FI        VARIABLE FORMAT FOR PUNCHED OUTPUT OF CONCENTRATION PROFILE.
C AMW,BMW   THE ATOMIC MASSES FOR THE PURE SPECIES OF ALPHA AND BETA.
C AMV,BMV   PARTIAL MOLAR VOLS. OF PURE SPECIES OF ALPHA AND BETA PHASES
C NDA,NDB,NDG NO. OF D VALUES IN THE ALPHA ,BETA AND GAMMA PHASES.
C CN        CONCENTRATION, ATOM FRACTION B.
C DN        DIFFUSION COEFFICIENT, IN CENTIMETER SQUARED PER SECOND.
C ND        NUMBER OF ( CONCENTRATION, DIFFUSION COEFFICIENT ) SETS.
C
C SPECIAL NOTE.. SOLUBILITY LIMITS FOR EACH PHASE ARE TAKEN
C AS THE END POINT CONCENTRATIONS (CN) FOR THAT PHASE.
C
C CPRIME    CONCENTRATION IN THE B RICH REGION.
C CZERO     CONCENTRATION IN THE MATRIX.
C           CPRIME,CZERO ARE USED TO SET UP A STARTING PROFILE
C RC        INITIAL LENGTH OF TOTAL DIFFUSION ZONE.
C RCF       FINAL LENGTH OF TOTAL DIFFUSION ZONE.
C RB        INITIAL THICKNESS OF BETA.
C RD        ESTIMATED DIFFUSION ZONE FOR BETA.
C RG        OUTER RADIUS OF INITIAL GAMMA PHASE.
C NXA,NXB,NXG NO. OF X-LOCATIONS IN THE ALPHA ,BETA AND GAMMA PHASES.
C
      READ (5,102)  ( ABLE(I), I=1,16 )
      READ(5,101) PCS
      READ (5,101)  NSP
      READ (5,101)  NO
      READ (5,102)(FI(I),I=1,8)
      READ (5,124)  AMW,BMW
      READ (5,124)    AMV,BMV
      READ (5,101)  NDA,NDB,NDG
      NDB1=NDB+1
      NDBG=NDB+NDG
      NDBG1=NDBG+1
      ND=NDB+NDG+NDA
      ND2=ND-2
      READ(5,110) (CN(I),DN(I),I=1,ND)
      READ (5,124)  CPRIME,CZERO
      READ(5,103) RC,RCF,RB,RD,RG
      READ (5,101)  NXA,NXB,NXG
C
C THIS SEGMENT INSPECTS THE TABLE OF DIFFUSION COEFFICIENTS TO SEE IF IT
C IS ARRANGED IN DESCENDING ORDER OF CONCENTRATION AND, IF NOT, ARRANGES
C IT SO.
C
      IF ( CN(1).GE.CN(ND) )  GO TO 4
      DO 2 I=1,ND
      CP(I)=CN(I)
    2 XSN(I) = DN(I)
      DO 3 I=1,ND
      N = ND+1-I
      CN(I) = CP(N)
    3 DN(I) = XSN(N)
    4 CONTINUE
C
C THE SOLUBILITY LIMITS ARE TAKEN AS FOLLOWS.
C
      CBG= CN(NDB)
      DBG= DN(NDB)
      CGB= CN(NDB1)
      DGB= DN(NDB1)
      CGA= CN(NDBG)
      DGA= DN(NDBG)
      CAG= CN(NDBG1)
      DAG= DN(NDBG1)
      DMAXA=0.
      DMAXB=0.
      DMAXG=0.
      DO 60 I=1,NDB
   60 IF(DMAXB.LT.DN(I)) DMAXB=DN(I)
      DO 61 I= NDB1,NDBG
   61 IF(DMAXG.LT.DN(I)) DMAXG=DN(I)
      DO 62 I=NDBG1, ND
   62 IF(DMAXA.LT.DN(I)) DMAXA=DN(I)
      WRITE(6,104)
      WRITE (6,105) ( ABLE(I), I=1,16 )
      WRITE(6,108) AMW,BMW,AMV,BMV
      WRITE(6,113) (CN(I),DN(I),I=1,ND)
      WRITE(6,118) CBG,DBG,CGB,DGB,CGA,DGA,CAG,DAG
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS IN ATOM FRACTION TO
C VOLUME FRACTION.
C
      DO 394 I=1,ND
      CN(I)=(CN(I)*BMV)/(AMV+CN(I)*(BMV-AMV))
  394 CONTINUE
      CZERO = CZERO *BMV /(AMV+CZERO*(BMV-AMV))
      CPRIME= CPRIME*BMV /(AMV+CPRIME*(BMV-AMV))
      CBG= CN(NDB)
      CGB= CN(NDB1)
      CGA= CN(NDBG)
      CAG= CN(NDBG1)
      CDELBG=CBG-CGB
      CDELGA=CGA-CAG
      T=0.
      IF (NSP.EQ.0)  GO TO 6
C
C THE FOLLOWING SEGMENT IS NOT USED UNLESS A STARTING PROFILE IS READ IN
C
C ABLE      IDENTIFICATION OF PROFILE.
C FI        VARIABLE FORMAT FOR READING IN PROFILE.
C T         TIME OF PROFILE IN SECONDS.
C ICAW      OPTION TO READ THE PROFILE IN ATOM FRACTION OR WT. FRACTION
C           +1, WEIGHT FRACTION
C           -1, ATOM FRACTION
C NP        NUMBER OF ( X, CONCENTRATION ) SETS IN PROFILE.
C CSN       CONCENTRATION, IN ATOM FRACTION B.
C XSN       DISTANCE, IN CENTIMETERS.
C
      READ (5,102)  ( ABLE(I), I=1,16 ),  ( FI(I), I=1,8 )
      READ (5,103)  T
      READ(5,101) ICAW
      READ (5,101)  NP
      READ (5,FI)   (CSN(I  ),XSN(I),I=1,NP)
      IF(ICAW.GT.0)GO TO 622
      DO 621 I=1,NP
  621   CP(I) =(CSN(I  )*BMV)/(AMV+CSN(I  )*(BMV-AMV))
      GO TO 623
  622 DO 620 I=1,NP
      CP(I)   =(CSN(I  )*(BMV/BMW))/((CSN(I  )*(BMV/BMW))+((1.0-CSN(I  )
     1)*(AMV/AMW)))
  620 CONTINUE
  623 CONTINUE
      CT=CDELBG/2.0
      ITEST=0
      KT=0
   29 KT=KT+1
      KT1=KT+1
      CTEST=CSN(KT  )-CSN(KT1  )
      IF (CTEST.LT.CT) GO TO 29
      IF(ITEST.EQ.1)GO TO 30
      RB =(XSN(KT)+XSN(KT1))/2.0
      ITEST=1
      IF(ITEST.EQ.1)CT=CDELGA
      GO TO 29
   30 RG=(XSN(KT)+XSN(KT1))/2.0
      WRITE (6,105)  ( ABLE(I), I=1,16 )
    6 TSTART=T
      TFINAL=T
      TSMIN=TSTART/60.0
      GEF= 4.* (1.+FLOAT(PCS))
C
C THE FOLLOWING SEGMENT SETS UP THE X ARRAY TO BE USED IN THE SOLUTION.
C NX IS THE NUMBER OF X-LOCATIONS WHERE THE COMPOSITION IS TO BE CALCU-
C LATED. NOTE THAT THERE ARE TWO POINTS AT EACH INTERFACE.
C
      NXB1=NXB+1
      NXB2=NXB+2
      NXB3=NXB+3
      NXBM1=NXB-1
      NXBM2=NXB-2
      NXBM3=NXB-3
      NXBG=NXB+NXG
      NXBG1=NXBG+1
      NXBG2=NXBG+2
      NXBG3=NXBG+3
      NXBGM1=NXBG-1
      NXBGM2=NXBG-2
      NX=NXB+NXG+NXA
      NXM1=NX-1
      THB=RB
      THG=RG-RB
      THA=RC-RG
      R1=RB-RD
      THBI=THB
      THGI=THG
      THAI=THA
      DXB=THB/(NXB-1)
      IF(R1.GT.0.0)DXB=(THB-R1)/(NXB-2)
      DXG=THG/(NXG-1)
      DXA=THA/(NXA-1)
      DMAX=AMAX1(DMAXA,DMAXB,DMAXG)
      TCONV=DMAX /RC**2
      XN(1) = 0.0
      XN(2)=DXB
      IF(R1.GT.0.0) XN(2)=R1
      DO 11  M=3,NXB
   11 XN(M) = XN(M-1) + DXB
      XN(NXB+1) =XN(NXB)
      DO 9 M=NXB2,NXBG
    9 XN(M)=XN(M-1)+DXG
      XN(NXBG1)=XN(NXBG)
      DO 72 M=NXBG2,NX
   72 XN(M)=XN(M-1)+DXA
      WRITE(6,114)NXB,NXG,NXA,THB,THG,THA,RC,RCF,RD,PCS
      IF(NSP.GT.0) GO TO 52
      DO 481 I=1,NXBM1
  481 CSN(I  )=CPRIME
      CSN(NXB  )=CBG
      CSN(NXB1  )=CGB
      DELC=(CGB-CGA)/(NXG-1)
      DO 39 I=NXB2,NXBG
   39 CSN(I  )=CSN(I-1  )-DELC
      CSN(NXBG1  )=CAG
      DO 41 I=NXBG2,NX
      CSN(I  )=CZERO
   41 CONTINUE
   52 DO 51 I=1,NX
   51 CP(I)= CSN(I)
      CONST=1.
      IF(PCS.EQ.1) CONST=3.14159
      IF(PCS.EQ.2) CONST= 4./3.*3.14159
      SUMI=0.
      DO 242 M=2,NX
  242 SUMI= SUMI+(CP(M)+CP(M-1)) *0.5*CONST*(XN(M) **(PCS+1)-XN(M-1)**
     1 (PCS+1))
      WRITE(6,104)
      IF( NSP.GT.0)  GO TO 411
      NGRID=0
      GO TO 410
C
C   CONTROL IS TRANSFERRED BACK TO THIS PT IF THE GRID REGION IS TO BE
C   EXPANDED IN EITHER THE BETA OR ALPHA PHASES
C
   12 DO 13 I=1,NX
   13 XSN(I)=XN(I)*RC
      NGRID=1
      WRITE(6,104)
      WRITE(6,107)TFINAL,THBF,THGF,THAF
      R1=R1*RC
      THG=THG*RC
      DXG=THG/(NXG-1)
      THB=THB*RC
      RD=THB-R1
      IF (CP(NX-2) .GT.CZERO1 ) RC=RC+RC-THB -THG
      IF( RC.GT.RCF) RC=RCF
      THA=RC-THB-THG
      IF(CP(4) .LT. CPRIME1) RD=RD*2
      R1=THB-RD
      IF(RD.GE.THB) R1=0.0
      DXB=(THB-R1)/(NXB-1)
      IF(R1.GT.0.0) DXB=(THB-R1)/(NXB-2)
      DXA=THA/(NXA-1)
      TCONV=DMAX /RC**2
      XN(1)=0.0
      XN(2)=DXB
      IF(R1.GT.0.0) XN(2)=R1
      DO 14 M=3,NXB
   14 XN(M)=XN(M-1)+DXB
      XN(NXB1)=XN(NXB)
      DO 10 M=NXB2,NXBG
   10 XN(M)=XN(M-1)+DXG
      XN(NXBG1)=XN(NXBG)
      DO 15 M=NXBG2,NX
   15 XN(M)=XN(M-1)+DXA
      IF( (THB/RC) .LE. 0.00001)  GO TO 53
  411 CSN(1  )=CP(1)
      DO 16 M=2,NXBM1
      I=1
   17 I=I+1
      IF(I.GT.NP) GO TO 22
      IF(XN(M).GT.XSN(I) ) GO TO 17
   22 IF (I.GT.NP) I=NP
      I1=I-1
   16 CSN(M)=CP(I1)-(CP(I1)-CP(I))/(XSN(I)-XSN(I1))*(XN(M)-XSN(I1))
      CSN(NXB)=CP(NXB)
      CSN(NXB1)=CP(NXB1)
   53 IF( (THG/RC) .LE. 0.00001)  GO TO 54
      DO 26 M=NXBG2,NXBGM1
      I=1
   27 I=I+1
      IF(I.GT.NP)GO TO 28
      IF(XN(M).GT.XSN(I))GO TO 27
   28 IF(I.GT.NP) I=NP
      I1=I-1
   26 CSN(M)=CP(I1)-(CP(I1)-CP(I))/(XSN(I)-XSN(I1))*(XN(M)-XSN(I1))
      CSN(NXBG)=CP(NXBG)
      CSN(NXBG1)=CP(NXBG1)
   54 IF( (THA/RC) .LE. 0.00001)  GO TO 410
      DO 18 M=NXBG2,NX
      IF(XN(M).GT.XSN(NX)) GO TO 20
      I=1
   19 I=I+1
      IF(XN(M).GT.XSN(I)) GO TO 19
      I1=I-1
      CSN(M)=CP(I1)-(CP(I1)-CP(I))/(XSN(I)-XSN(I1)) *(XN(M)-XSN(I1))
      GO TO 18
   20 CSN(M  )=0.0
   18 CONTINUE
  410 CONTINUE
      DO 21 M=1,NX
   21 CP(M)=CSN(M)
      IF(NGRID.EQ.0) WRITE(6,132)
      IF(NGRID.EQ.1) WRITE(6,133)
      WRITE(6,115)
      WRITE(6,111) (I,XN(I), CSN(I),I=1,NX)
      DO 24 M=1,NX
   24 XN(M)=XN(M)/RC
      E1OL=THB/RC
      E2OL=(THB+THG)/RC
      THB=THB/RC
      THG=THG/RC
      THA=THA/RC
      DXB=DXB/RC
      DXG=DXG/RC
      DXA=DXA/RC
      R1=R1/RC
      IF(NGRID.GT.0) GO TO 23
      NP=NX
      CPRIME1= CPRIME-0.0001
      CZERO1=CZERO+ 0.0001
C
C NSOL      NUMBER OF DIFFUSION TIMES AT WHICH THE OUTPUT IS REQUIRED.
C ATIME     DIFFUSION TIME IN SECONDS AT WHICH OUTPUT IS REQUIRED.
C           THE PROGRAM RECYCLES BACK TO STATEMENT 38 TO READ A NEW
C           DIFFUSION TIME. NSOL DIFFUSION TIMES ARE READ IN ALTOGETHER.
C
      READ(5,101) NSOL
      LSOL = 0
   38 LSOL = LSOL+1
      KZ=0
      KT=1
      READ(5,125) ATIME
   23 ITEST=0
C
C  FOLLOWING ARE SOME OF THE SPECIAL FEATURES OF THE PROGRAM SO AS TO
C  OPTIMISE THE COMPUTER TIME AND THE SOLUTION ACCURACY.
C  IF ANY OF THE PHASE THICKNESSES IS LESS THAN 1 PERCENT OF THE TOTAL
C  THICKNESS, DO NOT RECOMPUTE THE CONCENTRATIONS FOR THAT PHASE. THIS IS
C  DONE BY SETTING THE PARAMETERS NAAA,NBBB OR NGGG TO ZERO.
C  IF ANY PHASE THICKNESS IS LESS THAN 0.001 PERCENT OF THE TOTAL THICKNESS,
C  SET THAT PHASE THICKNESS TO ZERO.
C  INITIALLY, THE GAMMA PHASE CAN GROW OR SHRINK BY LESS THAN 2 PERCENT OF
C  ITS PREVIOUS THICKNESS.
C
   40 KZ=KZ+1
      CALL COEF
      TIA=0.25*(DXA**2)/DMAXA
      TIB=0.25*(DXB**2)/DMAXB
      TIG=0.25*DXG**2/DMAXG
      NAAA=1
      NGGG=1
      NBBB=1
      IF((CSN(NXBG1)-CSN(NX))  .LE.  1.E-06)  NAAA=0
      IF((CSN(NXB1)-CSN(NXBG))  .LE.  1.E-06)  NGGG=0
      IF((CSN(1)-CSN(NXB))  .LE.  1.E-06)  NBBB=0
      IF(E1OL.LE. .01) NBBB=0
      IF(THG.LE. .01) NGGG=0
      IF(E2OL.GE..99) NAAA=0
      IF(NAAA.EQ.0)  TIA= 1.E+70
      IF(NBBB.EQ.0)  TIB= 1.E+70
      IF(NGGG.EQ.0)  TIG= 1.E+70
      TI=AMIN1(TIA,TIB,TIG) *DMAX
      DE1OL=0.
      IF(E1OL.LE. .00001 .OR. E1OL .GE. .99999)  GO TO 68
      AGB= SQRT(DGB*DMAX*DP(NXB2   ))
      ABG= SQRT(DBG*DMAX*DP(NXBM1  ))
      DE1OL=1.0/CDELBG*(AGB/DMAX *(-CSN(NXB3  )+4.0*CSN(NXB2  )-3.0*CGB
     1)/2.0/DXG-ABG/DMAX *(CSN(NXBM2  )-4.0*CSN(NXBM1  )+3.0*CBG)/2.0/
     2DXB)
   68 DE2OL=0.
      IF(E2OL.LE. .00001 .OR. E2OL .GE. .99999)  GO TO 69
      AAG= SQRT(DAG*DMAX*DP(NXBG2  ))
      AGA= SQRT(DGA*DMAX*DP(NXBGM1 ))
      DE2OL=1.0/CDELGA*(AAG/DMAX *(-CSN(NXBG3  )+4.0*CSN(NXBG2  )-3.0
     1*CAG)/2.0/DXA-AGA/DMAX *(CSN(NXBGM2  )-4.0*CSN(NXBGM1  )+3.0*CGA)
     2/2.0/DXG)
   69 IF( NBBB.EQ.0 .OR. NAAA.EQ.0)  GO TO 70
      TEST= TI*(DE2OL-DE1OL)/THG
      IF(TEST.GT.0.02)TI=THG/(DE2OL-DE1OL)* 0.02
      IF(TEST.LT.-0.02)TI=-THG/(DE2OL-DE1OL)* 0.02
   70 CONTINUE
C
C   DE1OL - DISTANCE THE BETA/GAMMA INTERFACE MOVED PER UNIT TIME.
C   DE2OL - DISTANCE THE GAMMA/ALPHA INTERFACE MOVED PER UNIT TIME.
C
      IF(NBBB.EQ.0)  GO TO 47
      CP(1)=CSN(1  )+TI*GEF/2.   *DP(1)*(CSN(2  )-CSN(1  ))/DXB**2
C
C   DO LOOP 46 PERFORMS THE HEART OF THE CALCULATIONS FOR THE BETA PHASE
C
      DO 46 M=2,NXBM1
      C1=XN(M)/E1OL*DE1OL*(CSN(M+1  )-CSN(M-1  ))/2.0/DXB
      DJ= DP(M)*(ALOG(DP(M+1))- ALOG(DP(M-1)))
      C2=      (DP(M)*(PCS/XN(M)*(CSN(M+1  )-CSN(M-1  ))/2.0/DXB+
     1(CSN(M+1  )-2.0*CSN(M  )+CSN(M-1  ))/DXB**2)+ DJ/
     22.0/DXB*(CSN(M+1  )-CSN(M-1  ))/2.0/DXB)
      CP(M)=CSN(M  )+TI*(C1+C2)
      CSN(M-1  )=CP(M-1)
   46 IF(CP(M).GT.1.0) CP(M)=1.0
      IF(THG.LE. .00001 .AND. E2OL .GE. .99999)  CP(NXB)= CSN(NXB)+
     1TI*2. *DP(NXB) * (CSN(NXBM1)-CSN(NXB)) /DXB**2
      CSN(NXB )= CP(NXB )
      CSN(NXBM1  )=CP(NXBM1)
   47 CONTINUE
C
C   DO LOOP 42 IS THE HEART OF CALCULATIONS FOR THE GAMMA PHASE.
C
      IF(NGGG.EQ.0)  GO TO 48
      IF(E1OL.LE. .00001) CP(NXB1)= CSN(NXB1)+TI*GEF/2. *DP(NXB1)*
     1(CSN(NXB2)-CSN(NXB1)) /DXG**2
      IF(E2OL.GE. .99999) CP(NXBG)=CSN(NXBG)+TI *2.*DP(NXBG) *
     1(CSN(NXBGM1)- CSN(NXBG)) /DXG**2
      DO 42 M=NXB2,NXBGM1
      KPRIME=(XN(M)-E1OL)/(E2OL-E1OL)
      C1=(DE1OL+KPRIME*(DE2OL-DE1OL))*(CSN(M+1  )-CSN(M-1  ))/2.0/DXG
      DJ= DP(M)*(ALOG(DP(M+1))- ALOG(DP(M-1)))
      C2=      (DP(M)*(PCS/XN(M)*(CSN(M+1  )-CSN(M-1  ))/2.0/DXG+(CSN(
     1M+1  )-2.0*CSN(M  )+CSN(M-1  ))/DXG**2)+DJ*
     2(CSN(M+1  )-CSN(M-1  ))/4/DXG**2)
      CP(M)=CSN(M  )+TI*(C1+C2)
   42 CSN(M-1  )=CP(M-1)
      GO TO 45
   48 DO 63 M=NXB2,NXBGM1
      CP(M)= CP(NXB1) - (CP(NXB1)-CP(NXBG)) /FLOAT(NXG-1) *FLOAT(M-NXB1)
   63 CSN(M-1)= CP(M-1)
   45 CONTINUE
      CSN(NXBG)= CP(NXBG)
      CSN(NXBGM1  )=CP(NXBGM1)
C
C   DO LOOP 43 IS THE HEART OF CALCULATIONS FOR THE ALPHA PHASE.
C
      IF(NAAA.EQ.0)  GO TO 49
      IF(E1OL.LE. .00001 .AND. THG.LE. .00001)  CP(NXBG1) =CSN(NXBG1)+
     1TI*GEF/2. * DP(NXBG1) *(CSN(NXBG2)- CSN(NXBG1)) /DXA**2
      DO 43 M=NXBG2,NXM1
      IF(CP(M-1).EQ.0.0) GO TO 421
      C1=(1.0-XN(M))/(1.0-E2OL)*DE2OL*(CSN(M+1  )-CSN(M-1  ))/2.0/DXA
      DJ= DP(M)*(ALOG(DP(M+1))- ALOG(DP(M-1)))
      C2=DP(M)*(PCS/XN(M)*(CSN(M+1  )-CSN(M-1  ))/2.0/DXA+(CSN(M+1  )-
     12.0*CSN(M  )+CSN(M-1  ))/DXA**2)+         DJ      *(CSN(M+1  )-
     2CSN(M-1  ))/4.0/DXA**2
      CP(M)=CSN(M  )+TI*(C1+C2)
      IF(CP(M).LT.0.) CP(M)=0.
  421 IF(CP(M-1).LT.1.0E-05) CP(M)=0.0
   43 CSN(M-1  )=CP(M-1)
      CP(NX)=CSN(NX  )+TI*2.0*DP(NX)*(CSN(NX-1  )-CSN(NX  ))/DXA**2
      CSN(NXM1  )=CP(NXM1)
      CSN(NX  )=CP(NX)
   49 CONTINUE
      DBGIM=TI*DE1OL
      DGAIM=TI*DE2OL
      IF(E1OL.LE. 0.00001)  GO TO 67
      R1=R1+R1/THB*DBGIM
   67 CONTINUE
      E1OL=E1OL+DBGIM
      E2OL=E2OL+DGAIM
      TDBGIM=-THBI+E1OL*RC
      TDGAIM=-THBI-THGI+E2OL*RC
      DBGIM=DBGIM*RC
      DGAIM=DGAIM*RC
      IF( E1OL.LE. .00001) E1OL=0.
      IF( E2OL.LE. .00001) E2OL=0.
      IF( E1OL.GE. .99999) E1OL=1.
      IF( E2OL.GE. .99999) E2OL=1.
      THB=E1OL
      THG=E2OL-E1OL
      THA=1.0-THB-THG
      DXB=THB/NXBM1
      IF(R1.GT.0.0) DXB=(THB-R1)/NXBM2
      DXG=THG/(NXG-1)
      DXA=THA/(NXA-1)
      XN(1)=0.0
      XN(2)=DXB
      IF(R1.GT.0.0)XN(2)=R1
      DO 44 M=3,NX
      XX=DXB
      IF(M.GT.NXB) XX=DXG
      IF(M.GT.NXBG) XX=DXA
      XN(M)=XN(M-1)+XX
      IF(M.EQ.NXB1)XN(M)=XN(NXB)
   44 IF(M.EQ.NXBG1)XN(M)=XN(NXBG)
      THBF=THB*RC
      THGF=THG*RC
      THAF=THA*RC
      TI=TI/TCONV
      TFINAL=TFINAL+TI
C     IF(KZ.EQ.KT)WRITE(6,106) NBBB, NGGG, NAAA,KZ
C     IF(KZ.EQ.KT) WRITE(6,103) DE1OL,DE2OL,E1OL,E2OL,DXB,DXG,DXA,
C    1THBF,THGF,THAF,TIB,TIG,TIA,TI,TFINAL
      IF(KZ.EQ.KT) KT=KT+5000
      IF(R1.GT.0.0.AND.CP(4).LT.CPRIME1) GO TO 12
      IF(RC.LT.RCF.AND.CP(NX-2).GT.CZERO1) GO TO 12
      ITEST=ITEST+1
      IF (TFINAL.GE.ATIME) GO TO 430
C
C AT THIS POINT, THE PROGRAM CYCLES BACK TO STATEMENT NUMBER 40 TO
C PERFORM THE CALCULATIONS FOR THE NEXT TIME STEP.
C
      GO TO 40
  430 CONTINUE
      DO 241 M=1,NX
      XN(M)=XN(M)*RC
  241 CONTINUE
      SUMF=0.
      DO 250 M=2,NX
  250 SUMF= SUMF+(CP(M)+CP(M-1)) *0.5*CONST*(XN(M) **(PCS+1)-XN(M-1)**
     1 (PCS+1))
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS BACK TO ATOM
C FRACTION.
C
      DO 256 M=1,NX
  256 CP(M)=CP(M)*AMV/(BMV-CP(M)*(BMV-AMV))
      TFMIN=TFINAL/60.0
      CALL PLOTPF(CP,XN,NX,NA,XALM,XAL,YAL,DXN,DYN,XMIN,YMIN,XIN,YIN,
     1NCX,XCAP,NCY,YCAP,TYTLE,RESULT,LSOL,SPS)
      WRITE (6,104)
      WRITE (6,105)  ( ABLE(I), I=1,16 )
      WRITE(6,116)TSTART,TSMIN,     TFINAL,TFMIN,     THBI,THGI,THAI,
     1THBF,THGF,THAF,RC,            TDBGIM,TDGAIM,SUMI,SUMF
      WRITE(6,115)
      WRITE(6,111) (I,XN(I), CP (I),I=1,NX)
      IF(NO.NE.2) GO TO 58
      PUNCH 102,   (ABLE(I),I=1,16) , (FI(I),I=1,8 )
      PUNCH 103,  TSTART,TSMIN,     TFINAL,TFMIN,     THBI,THGI,THAI,
     1THBF,THGF,THAF,RC,            TDBGIM,TDGAIM,SUMI,SUMF
      PUNCH 101, NX
      PUNCH FI,   (CP(I),XN(I),I=1,NX)
   58 DO 260 I=1,NX
      XN(I)= XN(I) /RC
  260 CP(I)= CSN(I)
      IF ( LSOL.NE.NSOL ) GO TO 38
  992 CONTINUE
      CALL NFRAME
      STOP
      END
      SUBROUTINE PLOTPF(CP,XN,NP,NA,XALM,XAL,YAL,DXN,DYN,XMIN,YMIN,
     1XIN,YIN,NCX,XCAP,NCY,YCAP,TYTLE,RESULT,LSOL,SPS)
C
C THIS SUBROUTINE IS WRITTEN TO PLOT THE CALCULATED CONCENTRATION PROFILES
C ON A CALCOMP PLOTTER WITH THE CDC COMPUTER SYSTEM AT NASA LANGLEY RESEARCH
C CENTER, HAMPTON, VA. TO RUN ON OTHER COMPUTING SYSTEMS, CHANGES IN THE PLOT
C STATEMENTS MAY BE REQUIRED.
C
      DIMENSION CP(1),XN(1),XCAP(1),YCAP(1),TYTLE(1),RESULT(1)
      DO 1 I=1,NP
      CP(I)=CP(I)*100.0
    1 XN(I)=XN(I)*1.0E4
      IF( NA.GT.1 .AND. LSOL.GT.1 )  GO TO 100
      TXIN=XIN
      WH=XAL/80.
      IF (WH.LT.0.08) WH=0.08
      IF (WH.GT.0.25) WH=0.25
      DXL=XAL/DXN
      DYL=YAL/DYN
      PA=3.0
    3 NX=XN(NP)/TXIN+0.95
      DXM=NX
      XDN=DXN
      IF(SPS.EQ.0.) XDN= AMAX1(DXN,DXM)
      XBL=DXL*XDN
      IF(XBL.GT.XALM) TXIN=PA/2.0*XIN
      PA=PA+1.0
      IF(SPS. EQ.0. AND. XBL.GT.XALM) GO TO 3
      NX=DXN+1.1
      NY=DYN+1.1
      TS=0.0
      TF=1.25*WH
      WK=1.4*WH
      YAL=AMIN1(YAL,9.0)
C
C          DRAW X-AXIS AND Y-AXIS
C
      CALL CALPLT(2.0,2.5,-3)
      CALL CALPLT(XBL,0.0,2)
      CALL CALPLT(XBL,YAL,2)
      CALL CALPLT(0.0,YAL,2)
      CALL CALPLT(0.0,0.0,2)
      VYP=-2.*WK
      CX=NCX
      CSP=XBL/2.-(CX/2.+1.)*WK
C
C          PLACE NUMBERS ALONG X-AXIS
C
      DO 5 I=1,NX
      RI=I-1
      TP=RI*DXL
      TV=XMIN+RI*TXIN+0.1
      IF (TV.GE.100.) VXP=TP-1.5*WK
      IF (TV.LT.100..AND.TV.GE.0.2) VXP=TP-WK
      IF (TV.LT.0.2) VXP=TP-0.5*WK
      CALL CALPLT(TP,TS,3)
      CALL CALPLT(TP,TF,2)
      IF(TV.LE.0.001) GO TO 5
      CALL NUMBER (VXP,VYP,WK,TV,0.,-1)
    5 CONTINUE
C
C          LABEL X-AXIS
C
      CALL SCRIBE(CSP,-7.*WH,1.4*WH,0.05,XCAP,0.,NCX,9)
      CY=NCY
      CSP=YAL/2.0-(CY/2.0+1.0)*WK
C
C          PLACE NUMBERS ALONG Y-AXIS
C
      DO 10 I=1,NY
      RI=I-1
      TP=RI*DYL
      VYP=TP-0.5*WK
      TV=YMIN+RI*YIN+0.1
      IF (TV.GE.100.) VXP=-4.*WK
      IF (TV.LT.100..AND.TV.GE.0.2) VXP=-3.*WK
      IF (TV.LT.0.2) VXP=-2.*WK
      CALL CALPLT(TS,TP,3)
      CALL CALPLT(TF,TP,2)
      IF(TV.LE.0.0) GO TO 10
      CALL NUMBER (VXP,VYP,WK,TV,0.,-1)
   10 CONTINUE
      XAL2=XBL/4.0
      XAL3=XAL2+2.0
C
C          LABEL Y-AXIS
C
      CALL SCRIBE(-7.*WH,CSP,1.4*WH,0.05,YCAP,90.,NCY,9)
      CALL SCRIBE(XAL2,4.75,1.40*WH,0.05,RESULT,0.,20,9)
      CALL SCRIBE(XAL2,4.50,1.40*WH,0.05,TYTLE,0.0,80,9)
  100 CONTINUE
      DO 15 I=1,NP
      X =( XN(I)-XMIN)*XBL/(TXIN*DXN)
      Y=( CP(I  )-YMIN)*YAL/(YIN*DYN)
      IF (Y.LT.0.) Y=0.
      IF (Y.GT.YAL) Y=YAL
      IF(X.GT.XBL) X=XBL
      IF (X.LE.0..OR.Y.LE.0.) CALL CALPLT(X,Y,3)
      IF(X.GT.0.0.AND.Y.GT.0.0) CALL CALPLT(X,Y,2)
   15 CONTINUE
      DO 2 I=1,NP
      CP(I)=CP(I)/100.0
    2 XN(I)=XN(I)/1.0E4
      CALL CALPLT(0.0,0.0,3)
      IF(NA.EQ.1)  CALL NFRAME
      RETURN
      END
      SUBROUTINE COEF
C
C DIFFUSION COEFFICIENT CALCULATION BY LOGERTHIMIC INTERPOLATION.
C
      COMMON/STORE/NX,CN( 51),DN( 51),ND,CP(501),DP(501),DMAX
      DO 4 M=1,NX
      DO 2 I=2,ND
      IF(CP(M).LT.CN(I)) GO TO 2
      DP(M)=DN(I)*(DN(I-1)/DN(I))**((CP(M)-CN(I))/(CN(I-1)-CN(I)))
      GO TO 4
    2 CONTINUE
    4 DP(M)= DP(M)/DMAX
      RETURN
      END
